An Introduction to the samplr Package

The samplr package provides a number of functions for sampling from continuous distributions.

Rejection Sampling

Rejection sampling is a type of exact simulation method in numerical analysis and works for any distribution in \(\rm I\!R^m\) with a density. The idea of rejection sampling is to uniformly randomly sample the support and keep samples in the region under the graph of the density function.

Specifically, the steps of rejection sampling are:

  1. select uniformly at random a quantile candidate in the support of the function
  2. compute the probability density of the quantile candidate
  3. select a critical value via uniform sampling between 0 and the global maximum density
  4. retain the quantile candidate as a sample if the quantile density is less than the critical value; otherwise reject the candidate and start over.

Rejection sampling has two main assumptions:

  1. \(f(x) \le C\) (global maximum density)
  2. \(\exists\) a, b such that \(P(a \le X \le b) = 1\) (finite)

Function projectq2a

The projectq2a function generates random deviates from a continuous random variable. The user supplies the probability density function of the continuous random variable, and the function utilizing rejection sampling to generate the deviates.

Uniform

Here is an example of using projectq2a to sample from the standard uniform distribution.

The uniform distribution has probability density function:

\[\text{pdf = }\begin{cases}\frac{1}{b - a} \text{for } x \in [a, b]\\0 \text{ otherwise}\end{cases}\] We first generate a bunch of deviates and then plot their density in teal Then the true density is overlayed as a yellow line.

Beta

Here is an example of using projectq2a to sample from a beta distribution with both shape parameters equal to 2.

The beta distribution has probability density function:

\[\text{pdf = }\begin{cases}\frac{x^{\alpha - 1}(1 - x)^{\beta - 1}}{B(\alpha, \beta)} \text{ for }0 \le x \le 1\\0 \text{ otherwise}\end{cases}\]

\[\text{where }B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}\text{ with }\alpha, \beta \gt 0\]

We first generate a bunch of deviates and then plot their density in teal. Then the true density is overlayed as a yellow line.

Custom

Here is an example of using projectq2a to sample from a custom distribution.

In this case, the custom density function is that of a beta with shape parameters 2 and 5.

\[\text{pdf = }\begin{cases}30x(1 - x)^4 \text{ for }0 \le x \le 1\\0 \text{ otherwise}\end{cases}\]

We first generate a bunch of deviates and then plot their density in teal. Then the true density is overlayed as a yellow line.

Function projectq3a

The projectq3a function generates random deviates from a continuous 2D distribution defined on a square. The user supplies the probability density function, and the function utilizes rejection sampling to generate the deviate pairs.

2D Uniform on unit square

We first demo sampling from a uniform distribution defined on the unit square. In order to do so, we must first define a joint probability density function which takes x and y as arguments for the two variables.

\[\text{pdf = }\begin{cases}\frac{1}{(b - a)^2} \text{for } x, y \in [a, b]\\0 \text{ otherwise}\end{cases}\]

Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.

Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.

2D Beta on unit square

We first demo sampling from a beta distribution defined on the unit square. In order to do so, we must first define a joint probability density function which takes x and y as arguments for the two variables.

\[\text{pdf = }\begin{cases}\frac{x^{\alpha - 1}(1 - x)^{\beta - 1} + y^{\alpha - 1}(1 - y)^{\beta - 1}}{2B(\alpha, \beta)}\text{ for }0 \le x, y \le 1\\0\text{ otherwise}\end{cases}\]

\[\text{where }B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}, \text{ with }\alpha, \beta \gt 0\]

Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.

Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.

2D Custom on a square

We next demo sampling from a custom joint distribution defined on the square from -1 to +1. Again, we first define the joint pdf function which takes x and y as quantile arguments.

\[\text{pdf = }\begin{cases}\frac{3}{8}(x^2 + y^2) \text{for } x, y \in [-1, +1]\\0 \text{ otherwise}\end{cases}\]

Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.

Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.